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Abstract. We report on our numerical implementation of fully relativistic hydrodynamics coupled to Ein- 
stein's field equations in three spatial dimensions. We briefly review several steps in our code development, 
including our recasting of Einstein's equations and several tests which demonstrate its advantages for numer- 
ical integrations. We outline our implementation of relativistic hydrodynamics, and present numerical results 
for the evolution of both stable and unstable Oppenheimer-Volkov equilibrium stars, which represent a very 
promising first test of our code. 

INTRODUCTION 

The physics of compact objects is entering a particularly exciting phase. New instruments, including X-Ray and 
Gamma-Ray satellites and the new neutrino observatories SNO and Super-Kamiokande, can now yield unprecedented 
observations of neutron stars and black holes. Perhaps most excitingly, the new gravitational wave detectors LIGO, 
TAMA, GEO and VIRGO promise to open a gravitational wave window to the Universe and make gravitational 
wave astronomy a reality (see, e.g., [1]). 

Simultaneously, the availability of computational resources at modern supercomputers makes the simulation of 
realistic astrophysical scenarios involving relativistic compact objects feasable. Several groups, including two "Grand 
Challenge Alliances" [2], have launched efforts to construct numerical codes capable of solving Einstein's equations 
with or without matter sources in three dimensions, and simulating the merger of black hole or neutron star binaries 
(see, e.g., [3-9]). Numerical simulations will be necessary to predict the gravitational wave form from such processes 
to increase the likelihood of detections and, ultimately, to extract physical information from observations. Even 
though much progress has recently been made (see [9] for the most recent developments), significant obstacles still 
remain. 

In this contribution, we report on our systematic approach towards constructing such a numerical code. We first 
review our formulation of Einstein's equations and describe several tests, both for vacuum spacetimes and analytical 
matter sources, which demonstrate its advantages for numerical integrations. We then outline our implementation 
of hydrodynamics, and present promising numerical test results. We adopt geometrized units (G = c = 1) and the 
convention that Greek indices run from to 3, while Latin indices only run from 1 to 3. 



EVOLUTION OF THE GRAVITATIONAL FIELDS 

Most numerical implementations of Einstein's equations adopt a Cauchy formulation based on the 3+1 formulation 
by Arnowitt, Deser and Misner ( [10], see, e.g., [11] for an alternative characteristic formulation). However, a 
straightforward implementation of these "undressed" ADM equations tends to develop instabilities even for the 
evolution of small amplitude gravitational waves on a flat background (compare [12]). Following Shibata and 
Nakamura [13] and the spirit of many earlier one and two-dimensional codes [14], we have recently developed a 
modification of the ADM equations which proves to be much more suitable for numerical implementations [15]. 

More specifically, we modify the ADM equations in two ways. First, we split the spatial metric 7^ into a conformal 
factor exp((/)) and a conformally related metric jij according to 



7y = e '^Jij 



(1) 



and evolve (j) and separately. We choose (j) such that det(7ij) = 1, and similarly split the extrinsic curvature into 
its trace and its trace- free part. This split separates the "radiative" variables from the "non-radiative" variables in 
the spirit of the "York-Lichnerowicz" decomposition [16]. 

In the second stage we introduce "conformal connection functions" 

as independent functions (compare [17]). Some of the mixed second derivatives of in the conformal, spatial Ricci 
tensor Rij can now be written as first derivatives of the f". As a result, Rij becomes a manifestly elliptic operator 
on the metric 7^-. The analogous technique was used for the four-dimensional Ricci tensor Rap as early as in the 
1920's to make Einstein's equations manifestly hyperbolic [18]. For more details of our formulation, including an 
evolution equation for the r\ the reader is referred to [15] (see also the mathematical analysis in [19], and further 
numerical applications in [20]). 

In [15], we tested this form of Einstein's equations for small amplitude gravitational waves, and found that it 
performs far better than a similar implementation of the original ADM formulation. In particular, we found that wc 
could evolve such waves with harmonic slicing without encountering growing instabilities, whereas the original ADM 
code crashed after about 35 light-crossing times for the same initial data. In [21], we inserted analytic matter sources 
on the right hand side of Einstein's equations and evolved the fields in their presence. This approach allows us to 
study the numerical properties of the field evolution in the presence of highly relativistic matter sources without 
having to solve the equations of hydrodynamics: "hydro-without-hydro" . We inserted the Oppenheimer-Snyder 
solution for a relativistic, static star to test the long term stability of the evolution code, and the Oppcnhcimcr- 
Volkov solution for the collapse of a sphere of dust to a Schwarzschild black hole. These simulations focus on the 
highly relativistic, longitudinal fields, and complement our earlier tests involving dynamical transverse fields in [15]. 
With the code having passed these tests, we are now implementing both coUisionless matter, which will be described 
elsewhere, and hydrodynamics, as described below, to evolve the matter self-consistently with the fields. 

RELATIVISTIC HYDRODYNAMICS 

For a perfect fluid, the stress-energy tensor can be written 

T"'^ = (Po + Poe + P)u"u'^ + Pg'^l^, (3) 

where po is the rest mass density, e the specific internal density, P the pressure, ?i" the fluid four velocity, and gap 
the four dimensional spacetime metric. We construct constant entropy initial data with a polytropic equation of 
state 

P = Kpl, (4) 

where F = 1 -|- 1/n and n is the polytropic index, and where we assume the polytropic constant K to be unity 
without loss of generality. During the evolution, we adopt the gamma-law relation appropriate for adiabatic flow, 

P = (F - l)poe. (5) 

Following [9,22], we write the equation of continuity 

(PO«");a = (6) 



and the equations of motion 



T"^.Q = (7) 



in the form 



p*,i + (p,u*),» = 0, (8) 
e*,t + (e,w'),» = 0, (9) 
{p*Ui),t + {p*UiV^)j = -ae^^P^i - ap^uPa^i + p*ui(3'' ^ + 

P*UlUrn 



(^2f-<^,i + i (f 1,7'=- + f ^,7'')) • (10) 
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FIGURE 1. Mass versus central density for a n = 1 polytrope (left panel). We chose the two configurations marked by the 

filled circles {pc = 0.2 and pc — 0.4) as initial data for our dynamical simulations. The central density as a function of time 
for the evolution of the two initial data sets (right panel). As expected, the pc = 0.4 configuration soon collapses, whereas 
the Pc = 0.2 configuration remains stable for several dynamical timescales. 



Here we have defined the auxihary quantities 

= ae^^uVo, (11) 

e* = ae^^u°(poe)'/'^, (12) 

Ui = {l + Te)ui, (13) 

M" = (l + re)w", (14) 

v' = u'/vP = -P' + 7*^Uj/u°. (15) 

Similar equations have been used by many other groups (e.g. [23] and references therein). Wc integrate equations (8) - 
(10) with an artificial viscosity scheme suggested in [3] (see [8] for implementations of more elaborate shock capturing 
schemes). Given p*, e* and Ui on a new timelevel, can be found iteratively from the normaUzation relation 
u"Ua = — 1, which yields 

The matter sources for the right hand sides of Einstein's equations can then be constructed from these variables. 



NUMERICAL RESULTS 

As a first test of our implementation of hydrodynamics, we adopt the Oppenheimcr-Volkov solution describing 
equilibrium neutron stars in spherical symmetry as initial data and evolve these dynamically. Constructing a sequence 
of such Oppenheimer-Volkov solutions for increasing central rest mass densities pc, the ADM mass M of the star 
takes a maximum Mmax at p^"*' (see the left panel in Figure 1). For central densities smaller than p"'*, the star is 
in stable equilibrium, while for pc > p"'* it is unstable and will collapse to a black hole. 

Wc adopt a polytropic equation of state with F = 2 (n = 1), for which p"'* = 0.32 and Afmax = 0.164 when 
K = 1. We choose as initial data the two configurations marked with filled circles in Figure 1; a stable configuration 
with pc = 0.2, M = 0.157 and R/M = 5.5, and an unstable configuration with p^ = 0.4, M = 0.162 and R/M = 4.4. 

In the right panel of Figure 1, we show the central density as a function of time for the two initial configurations. We 
performed these runs on quite modest numerical grids with (64)"^ gridpoints using cartesian coordinates and imposing 
outgoing wave boundary conditions at x,y,z = 2. We employ harmonic slicing and zero shift. As expected, the 
unstable configuration soon collapses, while the stable configuration remains stable for several dynamical timescales 
(the period of the fundamental radial oscillation for a, Pc = 0.2 configuration is approximately 7pc , compare [9]). 
Ultimately, accumulation of numerical error causes this configuration to collapse too, but this can be delayed by 
increasing the grid resolution. 




In Figure 2 we show density contours at the beginning and towards the end of the simulation for the unstable 

configuration with pc = 0.4 initially. We also include arrows indicating the fluid flow Ui. The star is rapidly 
contracting and collapsing to a black hole. Up to the late stage of the collapse, the mass M is conserved to about 5 
%. We can follow the collapse to about a 18 fold increase of the central density. By this time, the central lapse has 
decreased from 0.4 initially to about 0.03. 



SUMMARY AND DISCUSSION 



We report on our systematic approach towards constructing a fully relativistic hydrodynamics code in three spatial 
dimensions. As part of this program, we have developed a new formulation of Einstein's equations which in several 
tests and applications has proved to be much more suitable for numerical implementations than the traditional ADM 
formulation [15,20]. Mathematical properties of this formulation have been analyzed in [19]. We have studied the 
evolution of small amplitude gravitational waves to test dynamical transverse fields and have inserted analytical 
solutions as matter sources (hydro-without-hydro) to test highly relativistic longitudinal fields. 

We outline our implementation of the relativistic equations of hydrodynamics and present preliminary test results 
for spherical neutron stars in hydrostatic equilibrium (Oppenheimer-Volkov stars). As expected, we find that stable 
stars remain stable for several dynamical timcscalcs, while unstable stars soon collapse to black holes. We conclude 
that our method seems like a very promising approach towards simulating the final plunge and merger of binary 
neutron stars. 

On a more speculative note, we point out that fully self-consistent hydrodynamics may not be a feasable approach 
towards simulating the coalescence and gravitational wave emission from neutron stars in the intermediate inspiral 
phase. In this epoch, the stars are very close and interact through a strongly relativistic tidal field, but reside 
outside the innermost stable circular orbit and hence move on a nearly circular orbit. Simulating the slow inspiral 
would require evolving the stars for hundreds or thousands of orbits, which is presently impossible. It may be 
possible, however, to insert known quasi-equilibrium binary configurations (e.g. [6,7]) into the field evolution code to 
get the transverse wave components approximately. Decreasing the orbital separation (and increasing the binding 
energy) at the rate consistent with the computed outflow of gravitational- wave energy would generate an approximate 
strong-field wave inspiral pattern. Such a "hydro-without-hydro" calculation may yield an approximate gravitational 
waveform from inspiraling neutron stars, without having to couple the matter and field integrations. 

Calculations were performed on SGI CRAY Origin2000 computer systems at the National Center for Supercom- 
puting Applications, University of Illinois at Urbana-Champaign. This work was supported by NSF Grants AST 
96-18524 and PHY 99-02833 and NASA Grant NAG 5-7152 at Illinois. 
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